Thermodynamics of the Antiferromagnetic Heisenberg Model on the Checkerboard 
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Employing numerical linked-cluster expansions (NLCEs) along with exact diagonalizations of fi- 
nite clusters with periodic boundary condition, we study the energy, specific heat, entropy, and 
various susceptibilities of the antiferromagnetic Heisenberg model on the checkerboard lattice. NL- 
CEs, combined with extrapolation techniques, allow us to access temperatures much lower than those 
accessible to exact diagonalization and other series expansions. We show that the high-temperature 
peak in specific heat decreases as the frustration increases, consistent with the large amount of un- 
quenched entropy in the region around maximum classical frustration, where the nearest-neighbor 
and next-nearest-neighbor exchange interactions (J and J', respectively) have the same strength, 
and with the formation of a second peak at lower temperatures. The staggered susceptibility shows 
a change of character when J' increases beyond 0.75J, implying the disappearance of the anti- 
ferromagnetic order at low temperatures. For J' = 4J, in the limit of weakly coupled crossed 
chains, we find large susceptibilities for stripe and Neel order with Q = (7r/2, it/2) at intermediate 
temperatures. Other magnetic and bond orderings, such as a plaquette valence-bond solid and a 
crossed-dimer order suggested by previous studies, are also investigated. 
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I. INTRODUCTION 



The checkerboard lattice is a unique two-dimensional 
(2D) system of great current interest. The next-nearest- 
neighbor (NNN) interactions, which are present on ev- 
ery other plaquette in a checkerboard pattern, not only 
can impose frustration and drive the system to exotic 
ground states but also provide a great tool for numer- 
ical and analytical investigators to study the evolution 
of physical properties in transitions between different ge- 
ometries. For instance, in the limit of weak NNN inter- 
actions, it is expected that the physics associated with 
the simple square lattice is dominant. In the antifer- 
romagnetic Heisenberg (AFH) model, this means a ten- 
dency toward long-range Neel ordering at temperatures 
smaller than the characteristic energy scale set by the 
nearest-neighbor (NN) magnetic exchange interaction, J. 
Whereas a ferromagnetic (negative) NNN exchange inter- 
action, J', favors this Neel ordering, an antiferromagnetic 
(positive) one introduces frustration and, thus, new types 
of ordering such as a valence-bond solid emerge. In the 
fully frustrated region where J ~ J' > 0, the lattice 
is a projection of the three-dimensional corner-sharing 
tetrahedrons (pyrochlore lattice) onto a 2D lattice. The 
other interesting limit is J' ^> J, where the 2D lattice 
is practically reduced to weakly coupled crossed chains, 
and physical properties are dominated by those of the 
one-dimensional (ID) system. Moreover, by eliminat- 
ing certain bonds, one can even turn the focus from the 
square basis of the underlying lattice to a triangular one 
that can capture the geometry of the Kagome lattice. 

The problem of the frustrated AFH model on the 
checkerboard lattice has its roots in early studies on 
its three-dimensional counterpart, the pyrochlore lattice. 
The latter system was originally studied by Harris et al^ 
using quantum field theory. They ruled out the possi- 
bility of a phase with long-range spin correlations but 



found strong correlation between NN spins, suggesting 
a dimerized ground state. A few years later, using per- 
turbative expansions and exact diagonalization, Canals 
and Lacroix? concluded that the ground state is a spin- 
liquid with correlations that decay exponentially by dis- 
tance. Around the same time, another study by Isoda 
and Mori)* in which a bond-operator approach was used, 
suggested a resonant-valence-bond-like plaquette phase. 

So far, the magnetic properties of the checkerboard lat- 
tice have been the focus of many theoretical studies* 4 - - — 
with compelling evidence that the ground state for J' = 
J (the planar pyrochlore) is a plaquette valence-bond 
solid (P-VBS) with long-range quadrumer order 4^ r 10 i 12 
This was shown by means of strong-coupling expan- 
sion,— £ exact diagonalization— as well as mean field 
theoryiS and a quadrumer boson approximation. 9 

In the limit of J' <C J, the existence of the long-range 
Neel order has also been established^^— Semiclassical 
approaches, such as the linear spin-wave ; 14 ' 15 and numer- 
ical resultsiii^ predict the stability of antiferromagnetic 
(AF) long-range order for J' /J < 0.75. However, this 
number is different in other studies that associate the in- 
stability of the P-VBS phase, as J' is reduced, with the 
transition to the Neel state (5/8 in Ref.[|, and 0.88-0.94 
in Ref.i). 

The situation in the limit of weakly coupled crossed 
chains ( J' 3> J) is less clear. There are at least two pro- 
posals for the ground state in this region of the parameter 
space; the first is the 2D spin-liquid ground state (sliding 
Luttinger liquid) characterized by the absence of long- 
range order and by elementary excitations being mass- 
less deconfined spinonsi 5 - This idea is supported by an 
exact diagonalization study of Sindzingre et a/., 13 which 
suggests a range of J/ J' = [0 — 0.8] for the ID behav- 
ior. However, their calculations suffer from strong finite- 
size effects even with 36 sites due to the quasi- ID nature 
of the problem. The second is the crossed-dimer (CD) 
phase suggested by Starykh et al£ They argued that, in 
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FIG. 1. Various ordered phases on the checkerboard lat- 
tice explored in this work; Neel order with a) Q = (tt,tv), 
b) Q = (7r/2,7r/2) (Neel*), c) Q = (n/2,ir), d) Q = (0,tt) 
(stripe). Open (solid) circles denote down-spins (up-spins); 
(e) crossed-dimer order where thick (thin) diagonal lines rep- 
resent strong (weak) bonds; and (f) P-VBS phase with strong 
dimer-dimer correlation between parallel bonds of uncrossed 
plaquettes marked by big circles. 



the CD phase, staggered dimer correlations, which have 
a power-law decay with distance in a perfect ID sys- 
tem, are stabilized when a weak interchain interaction 
(J) is present. As depicted in Fig. [He), in this phase, 
the "strong" (positive) dimers from perpendicular chains 
meet at the same crossed plaquette. This scenario is in 
agreement with the results of Arlego et fji.,— who exam- 
ined this idea by means of series expansion in terms of 
J and J' connecting the blocks of crossed dimers. In- 
cluding results from other works, Starykh et al£ also 
mapped out the global zero-temperature phase diagram 
of the system with respect to the ratio of J and J' and 
discussed the possibility of a magnetically ordered phase 
being present in the transition between the CD phase 
and the P-VBS phase. This so-called Neel* phase is the 
long-range ordered phase with diverging susceptibility at 
Q = (7t/2,7t/2) [see Fig. Ufa)]. Most recently, using a 
two-leg ladder to construct the 2D lattice in a density ma- 
trix renormalization group study, and by measuring var- 
ious spin-spin correlations, Moukouri 17 confirmed most 
of these predictions for the phase diagram except that 
the magnetically ordered phase in the proximity of the 
CD phase has a wave vector Q = (tt/2,tt) instead of the 
Q = (tt/2, 7t/2) proposed in Ref. [§■ Sketches of the for- 
mer order, along with the other orders explored here, are 
shown in Fig. [TJ 

Most of the numerical calculations for the AFH model 
on the checkerboard lattice have been done at zero tem- 
perature using finite clusters with periodic boundary con- 
dition. As discussed above, some of the early work o 7 ' 13 i 18 
helped shape theories that describe ground-state prop- 
erties such as the P-VBS. However, a systematic study 
of finite-temperature properties in the thermodynamic 
limit, more relevant to experiments, has been missing. 
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FIG. 2. Clusters generated in the first four orders of NLCE 
with a square building block on the checkerboard lattice. 



Our goal in this study is to explore the thermody- 
namic properties of this model and address the finite- 
temperature behavior of the susceptibilities to the or- 
dered phases proposed for the ground state and described 
above. 

We employ the numerical linked-cluster expansions 
(NLCEs) ) 20 i 21 along with exact diagonalization of finite 
clusters, to calculate thermodynamic properties of the 
system in different regions of the parameter space. We 
study the change in behavior of energy, entropy, spe- 
cific heat, and several susceptibilities as J and J' vary. 
We find that the high-temperature peak in specific heat 
is strongly suppressed in the case of maximum classical 
frustration, J' = J. Consistently, we see large amounts of 
unquenched entropy in this region, signaling the possibil- 
ity of a second peak in specific heat. Our study of differ- 
ent susceptibilities includes the staggered susceptibility, 
which for J'/ J < 0.75 continues to grow as the tempera- 
ture is lowered, suggesting that the ground state is Neel 
ordered with Q = (71, 7r) in this region. We also study 
the susceptibility to the the P-VBS phase using relevant 
order parameters and find that it is largest for J' ~ J. 
In the limit of weakly coupled crossed chains, and down 
to the lowest temperatures we can access, the dominant 
correlations belong to the Neel* and stripe phases. 

The paper is organized as follows: In Sec. [Til we present 
the model and briefly discuss NLCEs, and the extrap- 
olation techniques, along with the clusters utilized in 
the exact diagonalizations. The results are presented in 
Sec. IIIIl and a summary and conclusions are provided in 
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Sec. El 



Order Number of sites Number of clusters 

1 1 

1 4 1 

2 7 1 

3 10 2 

4 12 1 

4 13 4 

5 15 1 

5 16 10 

6 17 1 
6 18 7 
6 19 23 
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FIG. 3. (Color online) Periodic clusters on the checkerboard 
lattice used in our finite-size exact-diagonalization calcula- 
tions. The number inside each cluster represents its size. 



often gaining access to regions where most of the inter- 
esting phenomena take place. More details about these 
extrapolations can be found in the following subsection 
and references therein. 

Depending on the symmetry of the lattice and prop- 
erties of the model, the generation of clusters in NLCEs 
can be done using different building blocks. These in- 
clude the usual bond expansion, site expansion, triangle 
or square expansions, etc^i In this paper, we focus on 
the square expansion, which offers a particularly conve- 
nient approach in constructing the checkerboard lattice, 
i.e., by tiling it with crossed squares. In this picture, the 
first order in the expansion has a single crossed square, 
the second order has two crossed squares, and so on. The 
first four orders, including the zeroth order with a single 
site, are shown in Fig. [3J 

In the square expansion, the maximum number of sites 
of a cluster in the nth order is 3n + 1. Also, the number of 
topologically distinct clusters increases exponentially as 
the order increases. The number of clusters of each size, 
which need to be considered up to sixth order, is shown 
in Table HI Note that, out of 31 clusters in the sixth or- 
der, 23 have 19 sites, 7 have 18 sites, and 1 has 17 sites. 
Since the clusters have open boundaries, no translational 
symmetries can be used to block-diagonalize the Hamil- 
tonian matrix. This restricts the calculations to sixth 
or fewer orders, where, by using the conservation of the 
total spin in the z direction, we have to diagonalize ma- 
trices with linear size as large as (g 9 ) = 92, 378. This is 
nearly impossible using serial LAPACK subroutines on 
single-processor machines given memory restrictions and 
the time needed for such huge diagonalizations. There- 
fore, most of the calculations have been performed on 
parallel computers using SCALAPACK routines. 

Where possible, we compare results from NLCEs to 
those from exact diagonalization of finite clusters with 
periodic boundary conditions (ED) to build intuition 
about the finite-size effects that might have influenced 
results of previous studies. These clusters, with 16, 18, 
and 20 sites, are shown in Fig. [31 We use translational 
symmetries that are allowed on the checkerboard lattice 
and are not prohibited by the symmetries of the order 
parameters in the broken symmetry cases. The largest 
matrix we had to diagonalize in this case was for the 



II. MODEL AND NUMERICAL APPROACH 
The Hamiltonian 

The AFH Hamiltonian can be written as 

where S; is the spin-1/2 vector at site i, and (i, j) denotes 
bonds between NN sites i and j. denotes bonds 

between NNN sites i and j on every other square in a 
checkerboard pattern. 

Numerical Linked-Cluster Expansions 

NLCEs are linked-cluster expansion methods which al- 
low one to calculate the partition function and other ob- 
servables, per lattice site, in the thermodynamic limit at 
finite temperatures. The information for these quantities 
at a given temperature is built up by calculating contri- 
butions from all the clusters, up to a certain size, that 
can be embedded in the infinite lattice. Unlike high- 
temperature expansions (HTEs), each cluster is solved 
exactly using full diagonalization algorithms. Hence, NL- 
CEs have a region of convergence which extends beyond 
that of HTEs. Depending on the type of ordering that 
occurs in the system at low temperatures, NLCEs can re- 
main converged down to surprisingly low temperatures. 
Examples of these can be seen in the case of geometrically 
frustrated magnetic systems such as the Kagome lattice, 
where there is no long-range magnetic ordering^ - — As 
in other series expansion approaches, we use extrapola- 
tion techniques to perform the summation of existing or- 
ders to further decrease the temperature of convergence, 



TABLE I. Size and number of topologically distinct clusters 
up to the sixth order of the square expansion. 
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FIG. 4. (Color online) Energy per site vs temperature for the 
AFH model on the checkerboard lattice with NN and NNN ex- 
change interactions J = 0.50 and J' = 1.00, respectively. The 
thin dashed and dot-dashed lines represent the bare NLCE 
sums up to the fifth and sixth orders of the square expan- 
sion. The solid line shows the average of the last two terms 
in the Euler and Wynn extrapolations with the shaded (yel- 
low) area representing the "confidence limit" where all the 
extrapolations lie. The unit of energy is J'. 



sums, we take the average of the last two terms in the 
Euler and Wynn sums (solid line) . All these four extrap- 
olations lie in the shaded (yellow) region which can serve 
as the "confidence limit." We refer to this region around 
the average as the error bar, although it by no means 
represents statistical error bars. Below the temperature 
where bare NLCE sums diverge, the extrapolations' av- 
erage is not guaranteed to be the exact result in the ther- 
modynamic limit. However, along with the error bars, it 
serves as an estimate of the desired quantity. 



III. RESULTS AND DISCUSSIONS 

Here we study thermodynamic properties such as total 
energy, entropy, specific heat, and several magnetic sus- 
ceptibilities for a range of parameters, sweeping different 
regions of the phase diagram, from the simple square lat- 
tice without the NNN interaction to near the ID limit 
where NNN interactions dominate. For most of these 
quantities, we show results for J' = 0.00, 0.25, 0.50, 0.75, 
and 1.00 when J = 1.00 and J = 0.75, 0.50, and 0.25 
when J' = 1.00. The unit of energy is set to max( J, J'). 



20-site cluster, which had a linear dimension of 36,956. 



Extrapolations 

Measurements from all the clusters of every NLCE or- 
der are grouped together before summing different or- 
ders either regularly (bare sums) or by using Euler— or 
Wynr>2i sequence extrapolation algorithms. (For a de- 
tailed description of these algorithms see Ref. [2l|.) In the 
Euler sum, one can choose to have bare sums up to a 
particular order before using the Euler algorithm for the 
remaining orders. Here, we apply the Euler sum to the 
last four, three, two, and one terms. We find that the one 
with three Euler sums is generally the best (more phys- 
ically sensible). In the Wynn sum, we can have one or 
two cycles of improvement, each eliminating two terms, 
leaving us with four and two terms, respectively, out of 
the initial six. Because of the small number of terms in 
the Wynn sum, we find that using only one cycle yields 
a more reliable outcome. Hence, unless otherwise men- 
tioned, we show results throughout this paper for the 
Wynn sum with one cycle and Euler sum for three terms. 

The behavior of these extrapolations can be seen in 
Fig. |4l where we show, as an example, the energy per 
site (E) versus temperature for J = 0.50 and J' = 1.00. 
We also include the bare sums up to the fifth and sixth 
orders, which start diverging around T = 0.4J'. As ex- 
pected, the results from the Euler and Wynn sums show 
a less divergent behavior and extend the region of conver- 
gence to lower temperatures. To have a rough estimate 
for energy at temperatures not accessible by bare NLCE 



Energy, Entropy, Specific Heat and Bulk 
Susceptibility 

The specific heat per site (C) provides valuable infor- 
mation about the state of the system in different regions 
of the parameter space. In Fig. [5l we show this quan- 
tity after extrapolations of NLCE results for a range of 
values of J'/ J. For comparison, results from ED with 
18 and 20 sites are also shown. The highest peak ap- 
pears for the simple square lattice with no frustration 
[see Fig. [3](a)] . One can see that the bare NLCE results 
for fifth and sixth orders (dashed and dot-dashed lines, 
respectively) start deviating at T ~ 0.8, where the anti- 
ferromagnetic correlations presumably exceed the linear 
size of our biggest clusters. However, the average extrap- 
olation captures a peak around T = 0.6. More inter- 
estingly, both ED curves depart from the exact curve at 
a temperature greater than J and show almost no im- 
provement by increasing the cluster size from 16 to 20, 
with a position of the peak which is at slightly higher 
temperature. (The 16-site results are not shown.) 

As J'/ J increases to 0.5 [Figs. [5jb) and (c)], the peak 
in specific heat broadens, its maximum value decreases, 
and the temperature at which the latter is reached also 
decreases. Due to the increase in frustration, the AF 
correlations are suppressed and ED more accurately pre- 
dicts the location of the peak while still overestimating 
its value. For the same reason, the convergence in the 
bare NLCE sums is extended from T ~ 0.8 for J' = 0.0 
to T ~ 0.5 for J' = 0.5. Further increasing J' to 
0.75 [Fig. [5jd)] changes these features qualitatively by 
strongly suppressing the peak. In ED, the peak is pushed 
to lower temperatures (T ~ 0.3) and the agreement with 
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FIG. 5. (Color online) Specific heat vs temperature for various 
J and J': (a-d) J' < J and (e-h) J' > J. For comparison, 
results from ED with 18 and 20 sites are shown. The first 
peak is captured for all cases after extrapolation. The NLCE 
results are cut off roughly where the error bars exceed 0.1. 
The unit of energy is set to max(J, J'). 



exact NLCE results can be seen down to lower tempera- 
ture (T ~ 0.5) where the bare NLCE sums also diverge. 
These observations are consistent with results from pre- 
vious studies that find a transition at zero temperature 
from the magnetically ordered Neel phase to a disordered 
phase for J' > 0.75 

As expected, the minimum peak value is seen for the 
fully frustrated case of Fig. [5je) where J' = J [see also 
Figs. HJe) and intf)] ■ Although ED is in good agreement 
with NLCEs for T > 0.5, one can see significant finite- 
size effects at lower temperatures between the 18- and 
20-site clusters. The integral of C/T for the temperature 
range shown for the average extrapolation curve only re- 
covers about half of the entropy at infinite temperature, 
whereas 88% is recovered for the case of Fig. [SJa) with 
no frustration. At T = 0.3, the specific heat shows the 
tendency to develop a second peak. This tendency can 
be seen in both the NLCE and the ED results and, along 
with the fact that there is a huge amount of unquenched 
entropy already at T ~ 0.3 [see Figs. [BJc) and E^d)], 
strongly suggests that there is a second peak in specific 
heat at T < 0.3. 

As the value of J/ J' decreases from 1, the peak in spe- 
cific heat, shown in Figs. EJf )-Efh) , increases again and 
the 2D system starts to behave more and more like a ID 
chain. This can be inferred from the dramatic finite-size 
effects in ED. While the 20-site cluster can recover the 
NLCE results with relatively good accuracy, the results 




FIG. 6. (Color online) The evolution of (a, b) energy, (c, 
d) entropy, (e, f) specific heat, and (g, h) bulk susceptibility 
per site as J and J' change. These results are taken from 
the average extrapolations of NLCE and are cut off where 
the error bars reach 10% or less. Circles in (a), (e), and (g) 
are the data from a large-scale quantum Monte Carlo (QMC) 
simulation for J' = (Ref. l25l) . Circles in (c) are the result 
of a direct integration of C/T over temperature using the 
QMC results in (e), plus an additive constant to recover the 
infinite-temperature entropy, i.e., to account for the missing 
low-temperature tail of the specific heat. The statistical error 
bars for the QMC are smaller than the symbols and are not 
shown. 



for the 18-site cluster start deviating from NLCE at tem- 
peratures as high as 2.0. This can be understood from 
the fact that in the limit of decoupled crossed chains, 
,/ = 0, the 18-site cluster contains six decoupled peri- 
odic chains, each consisting of only three sites, whereas 
the 20-site cluster, contains two 10-site decoupled chains. 
Note that not only are the ID decoupled chains in the 
18-site cluster significantly smaller, but they also contain 
an odd number of chain sites, i.e., AF correlations are 
geometrically frustrated, and this strongly affects the re- 
sults. Due to the quantum fluctuations, any long-range 
order is suppressed near the ID limit, and so the extrapo- 
lations capture the specific heat with much smaller error 
bars for J' /J — 4 as seen in Fig. [3] (h). 

In Fig. [HI we show the evolution of energy, entropy (5), 
specific heat, and uniform susceptibility (x) per site as 
the value of J'/ J changes. One can see that the energy 
per site at temperatures below 0.2 increases monotoni- 
cally as J' increases and, as expected from the results in 
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tan _1 (J'/J)/jt 



FIG. 7. (Color online) Entropy divided by temperature as 
a function of frustration angle, tan -1 (J'/ J). A peak in en- 
tropy develops at J = J' as temperature is lowered below 
0.5. The values and the error bars are taken from the average 
extrapolation of NLCE results. 
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FIG. 8. (Color online) Log-log plot of extrapolated NLCE 
results for staggered susceptibility vs temperature. By intro- 
ducing J', the staggered susceptibility is suppressed. When 
J = J', X s ' 9 is more than one order of magnitude smaller 
than in the case of J' = 0. In the Euler extrapolation, only 
the last two terms have been used. The thin dotted lines are 
bare NLCE sums up to fifth and sixth orders. The thin solid 
lines represent the results from ED with 20 sites, and circles 
represent large-scale QMC results for J' = (Ref. [25). The 
statistical error bars for the QMC results are smaller than the 
symbols and are not shown. 



Order Parameter Susceptibilities 



Fig.0 the low-temperature entropy [Figs. EI c) andHJd)] 
is maximal in the case of J' = J. The previously dis- 
cussed decrease of the maximum value of the specific heat 
by increasing J' to J, followed by an increase for larger 
values of J' > J, is more clearly seen in Figs. Ete) and 
E^f). Finally, Figs. [6^g) andEJh) show that the uniform 
susceptibility remains small in all regions with a down- 
turn below T = 1.0. We have also included results from 
a large-scale stochastic series expansion QMC simulation 
(circles) with up to 256 x 256 spins for the unfrustrated 
case of J' = 25 using directed loop updatesi 26 ' 27 This is 
the only case that we consider where the low-temperature 
QMC calculation is not limited by the sign problem. 

To better compare the behavior of the entropy in dif- 
ferent regions, in Fig. [71 we show the entropy divided by 
temperature as a function of the frustration angle defined 
as (j> = tan _1 (J'/J). By lowering the temperature below 
0.5, the entropy develops a peak at J' = J which persists 
down to the lowest accessible temperature (with reason- 
able error bars for all angles). In the square lattice limit, 
<f> = 0, the specific heat, and therefore the entropy, are 
known to be quadratic in T at low temperatures. As can 
be seen in this figure, our results are consistent with this 
finding for T < 0.5. However, by increasing J'/ J to 1.0, 
this behavior changes completely and entropy decreases 
even more slowly than T. On the other hand, close to 
the ID limit, (f> > OAtt, the entropy has a linear region 
around T = 0.5 below which it decreases faster than T, 
similar to the weakly-frustrated regions with small <j). 



Other than the uniform susceptibility, which can be 
measured directly from the fluctuations of the total spin 
in the z direction, other susceptibilities per site are cal- 
culated using their definition as the second derivative of 
the free energy with respect to the field that couples to 
the corresponding order parameter (0): 



X 



T d 2 \nZ 



N dh 2 



h=0 



(2) 



where N is the number of sites, Z is the partition func- 
tion, and h is the field that couples to the order parameter 
in the new Hamiltonian, H' = H — hO. For example, we 
consider the following order parameter for Neel orderings 
with different wave vectors: 



O 



Nccl 



£V Q - R S*(R) 



(3) 



R 



where Q = (q x ,q y ), R runs over a Bravais lattice with 
the basis a = (^,0) and b = (0, f-), and S Z (R) is the 
total spin in the z direction in the corresponding unit 
cell. 

We find that the staggered susceptibility, x st9 [Q = 
(tt, 7t)], at low temperatures changes character when J' / J 
is increased from 0.75 to 1.00. As can be seen in Fig. [SI 
X st9 continues to grow by decreasing temperature in the 
weakly-frustrated region and as long as J'/ J < 0.75, but 
it shows a downturn at low T for J'/ J > 1.00. This is 
more clearly seen in the inset of Fig. [51 where we have 
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FIG. 9. (Color online) ED results for the susceptibility to 
the P-VBS order [see Eq. Q] per site vs temperature. Thick 
(thin) lines are results for the 20-site (16-site) cluster. The 
order is depicted in Fig. Q](f). The inset shows the suscepti- 
bility to the two-spin version of the plaquette order parameter 
as presented in Eq. ([5]). 
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FIG. 10. (Color online) Susceptibility to the crossed-dimer 
order [see Eq. (?? and Fig.QJe)] per site vs temperature. Thin 
dotted lines are the last two orders of bare NLCE sums, and 
thin solid lines are the ED results with 20 sites. In the Euler 
extrapolation, only the last two terms have been used. In the 
inset, lines are as in Fig[9]with thick (thin) lines representing 
ED results for the 20-site (16-site) cluster. 



plotted the inverse of the staggered susceptibility ver- 
sus temperature, and is consistent with previous find- 
mgg^Sr— which suggest that, in the latter region, the 
system no longer exhibits long-range Neel order. Note 
that the calculation of the staggered susceptibility for 
the unfrustrated case of J' = is one of the worst-case 
scenarios for NLCEs because the antiferromagnetic corre- 
lation length grows exponentially by decreasing the tem- 
perature. This can be realized by comparing the NLCE 
curve to the finite-size-converged (thermodynamic limit) 
QMC results (circles). Similar to the specific heat, the 
NLCEs results start deviating from the exact solution 
around T = 0.8. Nevertheless, NLCE provides a far bet- 
ter estimate for this quantity at low temperatures than 
ED. 

According to the Mermin- Wagner theorem,^? the 
Heisenberg model with finite-range exchange interactions 
in two dimensions cannot undergo a phase transition 
to a long-range ordered state at finite temperature by 
breaking a continuous symmetry. However, in light of 
the recent analytical and numerical predictions for the 
ground-state phases of this system, we calculate the 
finite-temperature susceptibilities associated with vari- 
ous order parameters to study their behavior as the tem- 
perature is lowered. These susceptibilities are shown in 
Figs. |9] to [12] in a low-temperature window for the rele- 
vant values of J and J' . 

In Fig. [9l we show the susceptibility to a plaquette or- 
der which is expected to be large in the P-VBS phase 
around J' = J. Fouet et alJ argued that the ground- 
state wave function in this phase is the symmetric com- 
bination of the pairs of singlets on parallel bonds of the 
uncrossed plaquettes. Based on that, we consider the 



following four-spin order parameter: 

4 = 32x^S a -S ;2 S, 3 .S Ml (4) 

O 

where I is the position of every other uncrossed square, 
marked by a circle in Fig.[T](f). The spin numbers around 
each of these squares are such that 1 and 2 (and there- 
fore, 3 and 4) are nearest neighbors. Since this kind 
of order involves uncrossed squares, NLCEs in crossed 
squares are not suited to measure the corresponding sus- 
ceptibility, and so we have obtained results only from 
ED. They show that this susceptibility is largest in the 
region around the maximum classical frustration. How- 
ever, significant finite-size effects are seen, especially for 
the J' = J case. In this region, the results for the 16-site 
cluster deviate from those for the 20-site cluster when 
T < 1.0, with the susceptibility being roughly a factor of 
2 larger at T ~ 0.1 for the 16-site cluster. Interestingly, 
for the 20-site cluster, the susceptibility shows a signif- 
icant decrease by further decreasing temperature below 
T = 0.07J. Note that most of the thermodynamic quan- 
tities, such as the specific heat and other susceptibilities 
calculated using ED (even with 20 sites), deviate from 
their exact NLCE counterparts (bare sums) starting from 
temperatures as high as 0.5 in this parameter region. So, 
the peak feature is expected to be a consequence of the 
finite-size nature of the calculations. We tested a more 
sophisticated order parameter suggested in Ref. to bet- 
ter capture the P-VBS phase, namely, the four-spin cyclic 
permutation operator (P4 + P^ 1 )^ and found the same 
qualitative results as for O4 after rescaling. 

Alternatively, one can define a simple two-spin order 
parameter as the sum of strong NN bonds around every 
other empty plaquette and weak NN bonds elsewhere to 
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FIG. 11. (Color online) Neel* susceptibility [Q = (7r/2,7r/2)] 
per site vs temperature for a range of ratios of J and J'. The 
inset shows the corresponding magnetic order where the open 
(solid) circles represent down-spins (up-spins) spins. Thin 
dotted lines are the last two orders of bare NLCE sums, and 
thin solid lines are the ED results with 16 sites. 

describe this phase: 

6 2 = X)(-!) i!e ( s ;i- s i2+S/2-S/3+S,3-S M +S;4-Su), (5) 

□ 

where I is the position of each empty square (□) in units 
of the NN lattice spacing and we have numbered the spins 
in each square clockwise, starting from the bottom left 
corner. The resulting susceptibilities for three values of 
J'/ J around the fully frustrated region are plotted in the 
inset of Fig. [9] and show similar trends as their four-spin 
counterparts. 

By decreasing J/J' to 0.5, we find that the low- 
temperature susceptibility to the CD order is enhanced 
(Fig. [ID)). To calculate this susceptibility, we consider the 
following order parameter: 

Ocd = 2 x ^(-r/*(S a • S, 3 + S /2 ■ S, 4 ), (6) 

SI 

where I is the position of each crossed square (EH) and 
spin numbering is the same as in Eq. ([5]) so that, Si 
and S3 (or S2 and S4) are at the two ends of diagonal 
bonds. Although this susceptibility is significantly larger 
for J' I J > 1, the extrapolated values for J'/ J = 4 ex- 
hibit a downturn at finite temperature. The results from 
ED with 20 sites overestimate the NLCE results at low 
T for J' > J. However, we see significant finite-size ef- 
fects between the 16- and 20-site clusters, shown in the 
inset of Fig. [TO] We have checked the susceptibility to 
a closely related order parameter in which there is one 
strong diagonal bond on every crossed plaquette, [specif- 
ically, Eq. ^ with a minus sign between the two terms] 
and found a behavior qualitatively similar the CD sus- 
ceptibility but with smaller values (not shown). Since 
the CD phase was predicted to exist for J' S> J 9 , an in- 
teresting question posed by these results is whether the 
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FIG. 12. (Color online) Susceptibility to the stripe order 
[Q = (0, 7r)] per site vs temperature. The inset shows the 
corresponding magnetic order where the open (solid) circles 
represent down-spins (up-spins) spins. Thin dotted lines are 
the last two orders of bare NLCE sums, and thin solid lines 
are the ED results with 20 sites. Circles are NLCE results for 
the Neel* order (Fig. [TT|) . 



peak feature will eventually disappear for smaller values 
of J I J' and one would find a susceptibility that always 
increases with decreasing temperature. In this scenario, 
the relevant temperature at which the CD phase becomes 
dominant is 0(J 2 / J') 9 , which is beyond the convergence 
region of our current NLCE calculations. 

We find that for large values of J' / J > 2 (weakly 
coupled crossed chains), there are two magnetic order- 
ings that are dominant at the lowest temperatures we 
can study. They are (i) the so-called Neel* order and 
(ii) stripes along the horizontal (or vertical) directions 
(Figs. [TT] and H"2"1) . The corresponding order parame- 
ters are defined in Eq. §5§ with Q = (7r/2,7r/2) and 
Q = (0,7r), and are depicted in Fig. QJb) andQJd), rc ~ 
spectively. The former has been proposed theoretically as 
the candidate for this region^ An intriguing observation 
is that the values for these two susceptibilities are hardly 
distinguishable, especially when J' > J. To illustrate the 
latter, we plot the NLCE results for the Neel* suscepti- 
bility against the stripe susceptibility in Fig. Q2] (circles) . 
One can see that the relative difference is negligible for 
all values of J'/ J shown. 

These results are consistent with what one would ex- 
pect at intermediate temperatures in the limit of J' 3> J, 
because both orders are compatible with the antiferro- 
magnetic correlations that develop along the diagonal 
chains. We note that for the ED with the 16-site cluster, 
using adjacency matrices, one can show that the mod- 
ified Hamiltonians, H' , are identical for the two order 
parameters. It would have been interesting to compare 
ED results for both orders with larger system sizes; how- 
ever, given the unit cell size for each order (eight sites for 
the Neel* and four sites for the stripe) and our computa- 
tional limitations with increasing system sizes, those re- 
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FIG. 13. (Color online) Susceptibility to the Neel order with 
Q = (-7T/2, 7r) per site, calculated using ED with 16 sites, 
versus temperature. The inset shows the corresponding order 
where empty (full) circles represent down (up) spins. 

suits are only available for the stripe order and are shown 
in Fig. 1121 Resolving which order becomes dominant at 
lower temperatures will require the study of larger cluster 
sizes, both in NLCEs and in ED. It is worth mentioning 
that, based on numerical calculations, the stripe order 
was suggested to be the one relevant to the ground state 
of the J i — J 2 model when J2 > 0.6Ji*i£ 

We have explored another magnetic ordering suggested 
by Mokouriii to be dominant in the limit J' ^> J. It has 
a wave vector of Q = (ir/2,%) and a unit cell of eight 
sites. Because of the breaking of certain symmetries of 
the lattice, the NLCE calculations for this case become 
much more expensive because one has to compute the 
physical properties of each cluster at different orienta- 
tions and locations on the lattice to properly deal with 
the broken symmetry. (This applies to the previously 
mentioned orders as well, but, due to the presence of 
other symmetries, computations are less costly in those 
cases.) Thus, we only present results from ED with 16 
sites for this type of order. As shown in Fig. Q21 not 
only is this susceptibility smaller close to the ID limit, 
but the maximum value, which belongs to the case of 
J' = </, is also much smaller than the maximum value 
seen for other orders with 16 sites (see, e.g., Fig. ITT1) . 
Therefore, a transition to this phase seems unlikely in 
any of the parameter regions. The fact that this type of 
order is not favored close to the ID limit is not surprising 
since, unlike in the Neel* or stripe ordered phases, spins 
on diagonal chains are not antiferromagnetically aligned. 



IV. SUMMARY 

We have calculated the thermodynamic properties of 
the AFH model on the checkerboard lattice using NL- 



CEs and ED, and studied their behavior as the system 
crosses over from a simple square lattice (J' = 0) to the 
maximally frustrated planar pyrochlore lattice ( J' = J) 
to the limit of one-dimensional crossed chains (J' ^S> J). 

We found that the peak value in the specific heat is sup- 
pressed as the frustration increases (by increasing J' / J 
from to 1), with strong indications that there is a sec- 
ond peak in the specific heat for J' ~ J. In the same 
region, finite-size effects in ED are minimal for tempera- 
tures above the convergence limit of NLCE. In contrast, 
close to the ID limit, ED results can vary significantly 
from one cluster to the other, depending on the size of 
periodic ID chains that exist inside each 2D cluster. Con- 
sistent with the reduced specific heat, entropy is maximal 
when J' = J at low temperatures with a decrease that is 
slower than T. 

We calculated the susceptibilities to several magnetic 
and bond orderings to explore the tendencies of the sys- 
tem toward different phases as the temperature is de- 
creased. By studying the staggered susceptibility, we 
found that the tendency toward Neel ordering with Q = 
(%, 7r) decreases appreciably when J'/ J > 0.75. By in- 
creasing the NNN interaction, antiferromagnetic correla- 
tions along the diagonal chains become important and 
other types of order emerge. To investigate this, we ex- 
amined the susceptibility to the P-VBS order using ED, 
and found that it is largest for J' ~ J. We also found 
large finite-size effects between 16- and 20-site clusters 
for J' = J. 

We further explored the susceptibility of the CD or- 
der, which is larger for J' > J but, according to the ex- 
trapolated NLCE results and for the values of J and J' 
considered here, does not dominate at the intermediate 
temperatures accessible within our NLCEs. Finite-size 
effects between the 16- and 20-site clusters were found to 
be significant in the ED calculations for J' > J. When 
J' > 2 J, i.e., for weakly coupled crossed chains, we found 
fast increasing susceptibilities at intermediate tempera- 
tures to Neel* order with Q = (tt/2, 7t/2), suggested by 
analytical results, and stripe order with Q = (0, it). Both 
of these orders are favored in this region due to the anti- 
ferromagnetically aligned spins along the chains. 
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